rm(list=ls())

library(foreign)
library(ggplot2)
library(dplyr)
library(stargazer)
library(sjPlot)
library(sjmisc)
library(gridExtra)
library(grid)
library(lattice)
library(ggpubr)
library(interflex)

options(scipen=999)



d_f1_f2 <- read.dta("final_observational_data.dta")

#Figure 1####
##Models for Figure 1####
#This is the output for Table A.1 in appendix
m1984 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
         data = d_f1_f2[d_f1_f2$year == "1984",])

m1988 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1988",])

m1992 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1992",])

m1996 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1996",])

m2004 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2004",])

m2008 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2008",])

m2012 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2012",])

m2016 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2016",])

m2020 <- lm(campact ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2020",])


##Preparing Plot for Figure 1####

#1984
#create the database for the first model
coef(m1984)
confint(m1984)
estimate <- 0.2004923184
ci_lo <- 0.0965624065
ci_hi <- 0.304422230
year <- 1984
plot_data <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))

#1988#
coef(m1988)
confint(m1988)
estimate <- 0.3050582187
ci_lo <- 0.0827305485
ci_hi <- 0.527385889
year <- 1988
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#1992
coef(m1992)
confint(m1992)
estimate <- 0.0533239381
ci_lo <- -0.0383644665
ci_hi <- 0.1450123428
year <- 1992
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#1996
coef(m1996)
confint(m1996)
estimate <- 0.1235687018
ci_lo <- 0.04100627287
ci_hi <- 0.206131131
year <- 1996
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2004
coef(m2004)
confint(m2004)
estimate <- 0.1050443416 
ci_lo <- -0.0330079983
ci_hi <- 0.243096682
year <- 2004
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2008
coef(m2008)
confint(m2008)
estimate <- 0.2161585730
ci_lo <- 0.11779567097
ci_hi <- 0.314521475
year <- 2008
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2012#
coef(m2012)
confint(m2012)
estimate <- 0.185691298
ci_lo <- 0.114089946
ci_hi <- 0.257292650
year <- 2012
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)


#2016
coef(m2016)
confint(m2016)
estimate <- -0.0514844598
ci_lo <- -0.168466799
ci_hi <- 0.065497879
year <- 2016
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2020
coef(m2020)
confint(m2020)
estimate <- 0.156731899
ci_lo <- 0.0671548843
ci_hi <- 0.246308913
year <- 2020
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

plot_data$year <- as.character(plot_data$year)
plot_data$JO <- c()

plot_data <- plot_data %>%
  group_by(year) %>% #do this analysis for each year, so want to group by year
  mutate(JO = case_when( 
    year == "1984" ~ "0", 
    year == "1988" ~ "0",
    year == "1992" ~ "1",
    year == "1996" ~ "1",
    year == "2004" ~ "1",
    year == "2008" ~ "0",
    year == "2012" ~ "0",
    year == "2016" ~ "1",
    year == "2020" ~ "0"
  ))


##Figure 1####
group.colors <- c("1" = "grey1", "0" = "grey0")


p <- ggplot(plot_data, aes(x = year, y = estimate, color = JO)) + #x is each treatment condition, y is the coefficient
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + #yintercept line at 0
  geom_point(aes(x = year,  #plotting the point
                 y = estimate, color = JO)) + 
  geom_linerange(aes(x = year, #plotting the 95% confidence interval
                     ymin = ci_lo,
                     ymax = ci_hi),
                 lwd = 1/2) + 
  theme_classic() +
  labs(y = "Campaign Activity",
       x = "Year") +
  scale_color_grey(start = 0.1, end = 0.6) +
  theme(legend.position="none")
p


#Figure 2####
##Preparing Data for Figure 2####
d_f1_f2 <- d_f1_f2 %>%
  na.omit() %>% #get rid of na, can't do the median if include na
  group_by(year) %>% #do this analysis for each year, so want to group by year
  mutate(hi_inc = case_when( #creae high income variable
    inc > median(inc) ~ 1, #if greater than the median then 1 for high income
    inc < median(inc) ~ 0 #if less than the median then 0 for high income
  )) %>%
  mutate(hi_inc_greater = case_when( #do the inverse of the above, should be exactly the opposite
    inc >= median(inc) ~ 1,
    inc < median(inc) ~ 0
  )) %>%
  mutate(hi_inc_less = case_when( #do the inverse of the above, should be exactly the opposite
    inc <= median(inc) ~ 1,
    inc > median(inc) ~ 0
  )) %>%
  ungroup() %>%
  as.data.frame()

##Models for Figure 2####
mf2_turnout <- glm(turnout ~ rgc*JO, data = d_f1_f2, 
                   family = binomial(link = "logit"))

summary(mf2_turnout)

mf2_donate <- glm(donate ~ rgc*JO, data = d_f1_f2, 
                  family = binomial(link = "logit"))

summary(mf2_donate)

##Table A.2####
stargazer(mf2_turnout, mf2_donate, type = "html", out = "Figure 2 Models.doc")

##Plot for Figure 2####
p2 <- plot_model(mf2_turnout, type = "pred", terms = c("rgc", "JO"),
                 title = "Voter Turnout (Low Cost)", colors = "bw") +
  theme_classic() +
  scale_y_continuous(labels = function(x) paste0(x/1), #i want a ratio not a percentage
                     limits = c(0, 1)) +
  ylab("Voter Turnout") +
  xlab("RGC") +
  scale_linetype(name = "",
                 labels = c("Non Black Candidate Years", "Black Candidate Years")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = 10))

p2

p3 <- plot_model(mf2_donate, type = "pred", terms = c("rgc", "JO"),
                 title = "Campaign Contributions (High Cost)", colors = "bw") +
  theme_classic() +
  scale_y_continuous(labels = function(x) paste0(x/1), #i want a ratio not a percentage
                     limits = c(0, 1)) +
  ylab("Campaign Contributions") +
  xlab("RGC") +
  scale_linetype(name = "",
                 labels = c("Non Black Candidate Years", "Black Candidate Years")) +
  theme(plot.title = element_text(hjust = 0.5)) +
  theme(plot.title = element_text(size = 10))

p3

ggarrange(p2, p3, ncol = 2, common.legend = TRUE, legend = "bottom")

#Figure 3####
##Models for Figure 3####
mf3a <- glm(donate ~ rgc*hi_inc_greater, data = d_f1_f2[d_f1_f2$JO == 0,], 
            family = binomial(link = "logit"))
summary(mf3a)

mf3b <- glm(donate ~ rgc*hi_inc_greater, data = d_f1_f2[d_f1_f2$JO == 1,], 
            family = binomial(link = "logit"))
summary(mf3b)

##Table A.3####
stargazer(mf3a, mf3b, type = "html", out = "Figure 3 Models.doc")


##Plot for Figure 3####


p4 <- plot_model(mf3a, type = "pred", terms = c("rgc [all]", "hi_inc_greater"),
                 title = "Non Black Candidate Years", colors = "bw") +
  theme_classic() +
  scale_y_continuous(labels = function(x) paste0(x/1), #i want a ratio not a percentage
                     limits = c(0, 0.75)) +
  ylab("Campaign Contributions") +
  xlab("RGC") +
  scale_linetype(name = "",
                 labels = c("Lower 50% Income", "Upper 50% Income")) +
  theme(plot.title = element_text(hjust = 0.5))

p4

p5 <- plot_model(mf3b, type = "pred", terms = c("rgc[all", "hi_inc_greater"),
                 title = "Black Candidate Years", colors = "bw") +
  theme_classic() +
  scale_y_continuous(labels = function(x) paste0(x/1), #i want a ratio not a percentage
                     limits = c(0, 0.75)) +
  ylab("Campaign Contributions") +
  xlab("RGC") +
  scale_linetype(name = "",
                 labels = c("Lower 50% Income", "Upper 50% Income")) +
  theme(plot.title = element_text(hjust = 0.5))

p5


ggarrange(p4, p5, ncol = 2, common.legend = TRUE, legend = "bottom")

#Figure A.1####
##These models make up table A.6 in appendix####
m1984s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1984",])

m1988s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1988",])

m1992s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1992",])

m1996s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "1996",])

m2004s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2004",])

m2008s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2008",])

m2012s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2012",])

m2016s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2016",])

m2020s <- lm(campactstandard ~ rgc + inc + female + polint + edu + age + relig, 
            data = d_f1_f2[d_f1_f2$year == "2020",])

#1984
coef(m1984s)
confint(m1984s)
estimate <- 0.767881040
ci_lo <- 0.369831831
ci_hi <- 1.165930249
year <- 1984
plot_data <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))

#1988
coef(m1988s)
confint(m1988s)
estimate <- 1.197118490
ci_lo <- 0.324653667
ci_hi <- 2.069583313
year <- 1988
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)
#1992
coef(m1992s)
confint(m1992s)
estimate <- 0.2047016762
ci_lo <- -0.147274768
ci_hi <- 0.556678120
year <- 1992
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#1996
coef(m1996s)
confint(m1996s)
estimate <- 0.668199298
ci_lo <- 0.22174193267
ci_hi <- 1.114656663
year <- 1996
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2004
coef(m2004s)
confint(m2004s)
estimate <- 0.702534959 
ci_lo <- -0.220756990
ci_hi <- 1.62582691
year <- 2004
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2008
coef(m2008s)
confint(m2008s)
estimate <- 1.172887610
ci_lo <- 0.6391654101
ci_hi <- 1.70660981
year <- 2008
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2012
coef(m2012s)
confint(m2012s)
estimate <- 0.94127983
ci_lo <- 0.578328476
ci_hi <- 1.30423118
year <- 2012
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)


#2016
coef(m2016s)
confint(m2016s)
estimate <- -0.3022657211
ci_lo <- -0.989070074
ci_hi <- 0.384538632
year <- 2016
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2020
coef(m2020s)
confint(m2020s)
estimate <- 0.890533567
ci_lo <- 0.381566734
ci_hi <- 1.3995004
year <- 2020
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#make year variable character
plot_data$year <- as.character(plot_data$year)
plot_data$JO <- c()

plot_data <- plot_data %>%
  group_by(year) %>% #do this analysis for each year, so want to group by year
  mutate(JO = case_when( 
    year == "1984" ~ "0", 
    year == "1988" ~ "0",
    year == "1992" ~ "1",
    year == "1996" ~ "1",
    year == "2004" ~ "1",
    year == "2008" ~ "0",
    year == "2012" ~ "0",
    year == "2016" ~ "1",
    year == "2020" ~ "0"
  ))

group.colors <- c("1" = "grey1", "0" = "grey0")


p_standardized <- ggplot(plot_data, aes(x = year, y = estimate, color = JO)) + #x is each treatment condition, y is the coefficient
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + #yintercept line at 0
  geom_point(aes(x = year,  #plotting the point
                 y = estimate, color = JO)) + 
  geom_linerange(aes(x = year, #plotting the 95% confidence interval
                     ymin = ci_lo,
                     ymax = ci_hi),
                 lwd = 1/2) + 
  theme_classic() +
  labs(y = "Campaign Activity",
       x = "Year") +
  scale_color_grey(start = 0.1, end = 0.6) +
  theme(legend.position="none")


p_standardized

#Figure A2####
##These models make up table A.7 in appendix####
m1984blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "1984",])

m1988blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "1988",])

m1992blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "1992",])

m1996blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "1996",])

m2004blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "2004",])

m2008blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "2008",])

m2012blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "2012",])

m2016blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "2016",])

m2020blf <- lm(campact ~ blf + inc + female + polint + edu + age + relig, 
             data = d_f1_f2[d_f1_f2$year == "2020",])

stargazer(m1984blf, m1988blf, m1992blf, m1996blf, m2004blf, m2008blf,
           type = "html", out = "black linked fate models.doc")

#1984
coef(m1984blf)
confint(m1984blf)
estimate <- 0.1046980421
ci_lo <- 0.0491878190
ci_hi <- 0.160208265
year <- 1984
plot_data <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))

#1988
coef(m1988blf)
confint(m1988blf)
estimate <- 0.1233265438
ci_lo <- 0.02891066
ci_hi <- 0.217742428
year <- 1988
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#1992
coef(m1992blf)
confint(m1992blf)
estimate <- 0.0712942153
ci_lo <- 0.03132580446
ci_hi <- 0.111262626
year <- 1992
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#1996
coef(m1996blf)
confint(m1996blf)
estimate <- 0.0348387553
ci_lo <- 0.00054920104
ci_hi <- 0.069128310
year <- 1996
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2004
coef(m2004blf)
confint(m2004blf)
estimate <- 0.0513308314 
ci_lo <- -0.016807186
ci_hi <- 0.119468849
year <- 2004
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2008
coef(m2008blf)
confint(m2008blf)
estimate <- 0.077775392
ci_lo <- 0.04003540852
ci_hi <- 0.115515376
year <- 2008
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2012
coef(m2012blf)
confint(m2012blf)
estimate <- 0.059956169
ci_lo <- 0.030744803
ci_hi <- 0.089167534
year <- 2012
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)


#2016
coef(m2016blf)
confint(m2016blf)
estimate <- 0.00349928197
ci_lo <- -0.057828735
ci_hi <- 0.064827299
year <- 2016
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#2020
coef(m2020blf)
confint(m2020blf)
estimate <- 0.068919179
ci_lo <- 0.0199710821
ci_hi <- 0.117867277
year <- 2020
plot_data2 <- as.data.frame(cbind(estimate, ci_lo, ci_hi, year))
plot_data <- rbind(plot_data, plot_data2)

#make year variable character
plot_data$year <- as.character(plot_data$year)
plot_data$JO <- c()

plot_data <- plot_data %>%
  group_by(year) %>% #do this analysis for each year, so want to group by year
  mutate(JO = case_when( 
    year == "1984" ~ "0", 
    year == "1988" ~ "0",
    year == "1992" ~ "1",
    year == "1996" ~ "1",
    year == "2004" ~ "1",
    year == "2008" ~ "0",
    year == "2012" ~ "0",
    year == "2016" ~ "1",
    year == "2020" ~ "0"
  ))

group.colors <- c("1" = "grey1", "0" = "grey0")


p_lf <- ggplot(plot_data, aes(x = year, y = estimate, color = JO)) + #x is each treatment condition, y is the coefficient
  geom_hline(yintercept = 0, colour = gray(1/2), lty = 2) + #yintercept line at 0
  geom_point(aes(x = year,  #plotting the point
                 y = estimate, color = JO)) + 
  geom_linerange(aes(x = year, #plotting the 95% confidence interval
                     ymin = ci_lo,
                     ymax = ci_hi),
                 lwd = 1/2) + 
  theme_classic() +
  labs(y = "Campaign Activity",
       x = "Year") +
  scale_color_grey(start = 0.1, end = 0.6) +
  theme(legend.position="none")

p_lf


#Table A.8####
#testing for multicollinearity, variance inflation factors
car::vif(m1984)
car::vif(m1988)
car::vif(m1992)
car::vif(m1996)
car::vif(m2004)
car::vif(m2008)
car::vif(m2012)
car::vif(m2016)
car::vif(m2020)



#Figure A3####
A17a <-ggplot(d_f1_f2, aes(x = rgc)) +
  geom_density(aes(group = as.factor(hi_inc_greater),
                   fill = as.factor(hi_inc_greater)), alpha = 1.2) +
  ylim(0, 3)+
  scale_fill_manual(name = "Income", 
                    labels = c("Low Income", "High Income"),
                    values = c("grey", "black"))+
  xlab("RGC")+  
  ylab("Density") +
  theme_classic() +
  theme(legend.position='bottom', legend.text=element_text(size=10),
        legend.direction = "horizontal") + 
  ggtitle("Income") +
  theme(plot.title = element_text(hjust = 0.5)) 
A17a

A17b <-ggplot(d_f1_f2, aes(x = rgc)) +
  geom_density(aes(group = as.factor(JO) , fill = as.factor(JO)), alpha = 1.2) + 
  ylim(0, 3)+
  scale_fill_manual(name = "Racial Context", 
                    labels = c("Non Black Candidate Year", "Black Candidate Year"),
                    values = c("grey", "black"))+
  xlab("RGC")+ 
  ylab("Density") +
  theme_classic() +
  theme(legend.position='bottom',  legend.direction = "horizontal",
        legend.text=element_text(size=10)) + 
  theme(plot.title = element_text(hjust = 0.5))  +
  ggtitle("Racial Context") 
A17b

densityplot <- ggarrange(A17a, A17b, ncol = 2)
densityplot


#Figure A4####
s6.linear <- interflex(estimator="kernel", method="linear", data = d_f1_f2, 
                       Y = "turnout", D = "hi_inc_greater", X = "rgc")
s6.plot <- plot(s6.linear)
s6.predict <- predict(s6.linear, pool = TRUE, 
                      subtitles=c('Low Income','High Income'),
                      color = c("grey", "black"))
turnoutplot<-s6.predict+ggtitle("Turnout")+
  theme(legend.position='bottom', legend.text=element_text(size=12))  +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) 

turnoutplot

s6.linear <- interflex(estimator='kernel', method='linear', data = d_f1_f2, 
                       Y = "donate", D = "hi_inc_greater", X = "rgc")
s6.plot <- plot(s6.linear)
s6.predict <- predict(s6.linear, pool = TRUE, 
                      subtitles=c('Low Income','High Income'),
                      color = c("grey", "black"))
donationplot<-s6.predict+ 
  ggtitle("Donations")+
  theme(legend.position='bottom', legend.text=element_text(size=10))+
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) 
donationplot

incomeplot <- ggarrange(turnoutplot, donationplot, ncol = 2, common.legend = TRUE, legend = "bottom")
incomeplot

#Figure A5####
s6.linear <- interflex(estimator='kernel', method='linear', data = d_f1_f2, 
                       Y = "donate", D = "JO", X = "rgc")
s6.plot <- plot(s6.linear)
s6.predict <- predict(s6.linear, pool = TRUE, 
                      subtitles=c('Non Black Candidate Years','Black Candidate Years'),
                      color = c("grey", "black"))
donationplotjo<-s6.predict+
  ggtitle("Donations")+
  theme(legend.position='bottom', legend.text=element_text(size=12))+
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5))
donationplotjo

s6.linear <- interflex(estimator='kernel', method='linear', data = d_f1_f2, 
                       Y = "turnout", D = "JO", X = "rgc")
s6.plot <- plot(s6.linear)
s6.predict <- predict(s6.linear, pool = TRUE, 
                      subtitles=c('Non Black Candidate Years','Black Candidate Years'),
                      color = c("grey", "black"))
turnoutplotjo<-s6.predict+
  ggtitle("Turnout")+
  theme(legend.position='bottom', legend.text=element_text(size=12)) +
  theme_classic() +
  theme(plot.title = element_text(hjust = 0.5)) 
turnoutplotjo

joplot <- ggarrange(donationplotjo, turnoutplotjo, ncol = 2, common.legend = TRUE, legend = "bottom")
joplot

#Figure A6####
out5 <- interflex(Y = "turnout", D = "hi_inc_greater", X = "rgc", Z = NULL, 
                  data = d_f1_f2, estimator = "binning", vcov.type = "robust", 
                  main = "Marginal Effects", ylim = c(0, 1), Xdistr = "density",
                  theme.bw = TRUE)
incvotebin<-plot(out5,  
                 ylab = "Marginal Effect of Income on Voting", 
                 main = "Voting DV")
incvotebin
out6 <- interflex(Y = "donate", D = "hi_inc_greater", X = "rgc", Z = NULL,
                  data = d_f1_f2, estimator = "binning", vcov.type = "robust",
                  main = "Marginal Effects", ylim = c(0, 1), Xdistr = "density",
                  theme.bw = TRUE)
incdonatebin<-plot(out6,  
                   ylab = "Marginal Effect of Income on Donations", 
                   main = "Donating DV")
incdonatebin
c<-grid.arrange(incvotebin, incdonatebin, nrow = 1)
c

#Figure A7####
out5 <- interflex(Y = "turnout", D = "JO", X = "rgc", Z = NULL, 
                  data = d_f1_f2, estimator = "binning", vcov.type = "robust", 
                  main = "Marginal Effects", ylim = c(0, 1), Xdistr = "density")
incvotebin<-plot(out5,  
                 ylab = "Marginal Effect of Black Candidate on Voting", 
                 main = "Voting DV")
incvotebin
out6 <- interflex(Y = "donate", D = "JO", X = "rgc", Z = NULL, 
                  data = d_f1_f2, estimator = "binning", vcov.type = "robust", 
                  main = "Marginal Effects", ylim = c(0, 1), Xdistr = "density")
incdonatebin<-plot(out6,  
                   ylab = "Marginal Effect of Black Candidate on Donations", 
                   main = "Donating DV")
incdonatebin
c<-grid.arrange(incvotebin, incdonatebin, nrow = 1)

library(e1071)
kurtosis(d_f1_f2$rgc)
library(lmom)
lmrd(samlmu(d_f1_f2$rgc))
ggsave(filename="kurtosis.jpeg", device="jpeg", height=5, width=10, units="in", dpi=500)
library(Lmoments)
Lmoments(d_f1_f2$rgc, rmax = 4, na.rm = FALSE, returnobject = FALSE, trim = c(0, 0))

#Table A.9#### 

d_f1_f2$rgclow <- 0
d_f1_f2$rgclow[d_f1_f2$rgc <= 0.688] <- 1
d_f1_f2$rgchigh <- 0
d_f1_f2$rgchigh[d_f1_f2$rgc >= 0.833] <- 1


difftest_lm <- function(x1, x2, model){
  diffest <- summary(model)$coef[x1,"Estimate"]-summary(model)$coef[x2,"Estimate"]
  vardiff <- (summary(model)$coef[x1,"Std. Error"]^2 + 
                summary(model)$coef[x2,"Std. Error"]^2) - (2*(vcov(model)[x1, x2])) 
  # variance of x1 + variance of x2 - 2*covariance of x1 and x2
  diffse <- sqrt(vardiff)
  tdiff <- (diffest)/(diffse)
  ptdiff <- 2*(1-pt(abs(tdiff), model$df, lower.tail=T))
  upr <- diffest + qt(.975, df = model$df)*diffse # will usually be very close to 1.96
  lwr <- diffest + qt(.025, df = model$df)*diffse
  df <- model$df
  return(list(est=round(diffest, digits =2), 
              t=round(tdiff, digits = 2), 
              p=round(ptdiff, digits = 4), 
              lwr=round(lwr, digits = 2), 
              upr=round(upr, digits = 2),
              df = df))
}

#diff tests
#donating DV, Black candidate
m1<-lm(donate ~ rgchigh*JO + rgclow*JO, data = d_f1_f2)
#voting DV, Black candidate
m2<-lm(turnout ~ rgchigh*JO + rgclow*JO, data = d_f1_f2)
#donating DV, income
m3 <- lm(donate ~ rgchigh*hi_inc_greater + rgclow*hi_inc_greater, data = d_f1_f2)
#voting DV, income
m4 <- lm(turnout ~ rgchigh*hi_inc_greater + rgclow*hi_inc_greater, data = d_f1_f2)

difftest_lm("rgchigh:JO", "JO:rgclow", m1)
difftest_lm("rgchigh:JO", "JO:rgclow", m2)
difftest_lm("rgchigh:hi_inc_greater", "hi_inc_greater:rgclow", m3)
difftest_lm("rgchigh:hi_inc_greater", "hi_inc_greater:rgclow", m4)

#Table A.10####
fit <- lm(donate ~ JO * hi_inc_greater * rgc, data = d_f1_f2)
summary(fit)
stargazer(fit, type = "html", out="tripleinteraction.doc")
